Source Detection and Measurement¶

Image
Contact authors: Douglas Tucker, Alex Drlica-Wagner
Last verified to run: 2024-06-02
LSST Science Pipelines version: Weekly 2024_16
Container Size: medium
Targeted learning level: Intermediate

Description: Access, display, and manipulate images; detect, deblend, and measure sources; and extract, plot, and use object footprints.

Skills: Run source detection, deblending, and measurement. Use source footprints.

LSST Data Products: DP0.2 processed visit images and catalogs.

Packages: lsst.pipe.tasks.characterizeImage, lsst.meas.algorithms.detection.SourceDetectionTask, lsst.meas.deblender.SourceDeblendTask, lsst.meas.base.SingleFrameMeasurementTask, lsst.meas.base.ForcedMeasurementTask

Credit: Originally developed by Alex Drlica-Wagner and Imran Hasan in the context of the LSST Stack Club. Forced photometry routine by Tobias Géron and Melissa Graham

Get Support: Find DP0-related documentation and resources at dp0.lsst.io. Questions are welcome as new topics in the Support - Data Preview 0 Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook provides a brief introduction to running the LSST Science Pipelines source detection, deblending, and measurement tasks. It does not go into depth about optimizing these tasks for different types of sources. It also covers how to access the footprint of a detected source and how to do forced photometry at user-specified coordinates.

Some source detection and measurement details come from Robert Lupton's Tune Detection.ipynb and Kron.ipynb. Interaction with lsst.afw.display was also improved by studying Michael Wood-Vasey's DC2_Postage Stamps.ipynb. More information on footprints can be found on the Stack Club notebook by Imran Hasan here.

1.1. Import packages¶

The numpy and matplotlib are widely used Python plot libraries.

From the LSST Science Pipelines, modules for data access via the butler, image display, image geometry, table handling, and various tasks for image characterization and source detection, deblending, and measurement are imported.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import getpass

import lsst.daf.base as dafBase
from lsst.daf.butler import Butler
import lsst.afw.image as afwImage
import lsst.afw.display as afwDisplay
import lsst.afw.table as afwTable
import lsst.geom as geom

from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask
from lsst.meas.algorithms.detection import SourceDetectionTask
from lsst.meas.deblender import SourceDeblendTask
from lsst.meas.base import SingleFrameMeasurementTask
from lsst.meas.base import ForcedMeasurementTask

1.2. Define functions and parameters¶

Set afwDisplay to use matplotlib as its backend. Set the matplotlib plot color table to be colorblind-friendly. Set the default figure size to be a square that will fit nicely into the notebook.

In [2]:
afwDisplay.setDefaultBackend('matplotlib')
plt.style.use('tableau-colorblind10')
plt.rcParams['figure.figsize'] = (8.0, 8.0)

Instantiate butler access to the DP0.2 data sets.

In [3]:
butler = Butler('dp02', collections='2.2i/runs/DP0.2')

2. Image access, display, and manipulation¶

2.1. Retrieve and display a calexp¶

Define the dataId for pre-selected visit 512055 and detector 75, then use the butler to retrieve the calexp (processed visit image).

In [4]:
dataId = {'visit': 512055, 'detector': 75}
In [5]:
calexp = butler.get('calexp', dataId=dataId)

A calexp contains the image plus a mask and variance planes.

Option: uncomment any of the following lines and execute the following cell to learn more about the various aspects of the calexp.

In [6]:
# calexp.maskedImage
# calexp.maskedImage.image
# calexp.maskedImage.mask
# calexp.maskedImage.variance
# calexp.image
# calexp.mask
# calexp.variance

The retrieved calexp also contains the point-spread function (PSF), the world coordinate system (WCS), and the photometric calibration for the image.

These components could be retreived with the following statements.

psf = calexp.getPsf()
wcs = calexp.getWcs()
pcal = calexp.getPhotoCalib()

Later steps in this tutorial perform source detection and measurement, so clear the existing DETECTED mask plane of the image.

In [7]:
calexp.mask.removeAndClearMaskPlane('DETECTED')

Display the calexp image.

In [8]:
fig = plt.figure()
afw_display = afwDisplay.Display()
afw_display.scale('asinh', 'zscale')
afw_display.mtv(calexp.image)

<Figure size 800x800 with 0 Axes>
Image

Figure 1: Above, the calexp image for visit 512055 and detector 75. Pixel count values range from -300 to 400.

2.2. Add the sky background back into the image¶

As an example of how to manipulate an image before reprocessing, use the butler to retrieve the sky background which was subtracted, and add it back in to the image.

In [9]:
bkgd = butler.get('calexpBackground', dataId=dataId)

Display the background.

In [10]:
fig = plt.figure()
afw_display = afwDisplay.Display()
afw_display.scale('linear', 'zscale')
afw_display.mtv(bkgd.getImage())
plt.title("Local Polynomial Background")
plt.show()

<Figure size 800x800 with 0 Axes>
Image

Figure 2: Above, the subtracted sky background for visit 512055 and detector 75. Pixel count values range from 3816 to 3822.

Add the background into the calexp and display the results.

Warning! Only execute this cell once or the background will be re-added multiple times.

In [11]:
calexp.maskedImage += bkgd.getImage()
In [12]:
plt.figure()
afw_display = afwDisplay.Display()
afw_display.scale('asinh', 'zscale')
afw_display.mtv(calexp.image)
plt.show()

<Figure size 800x800 with 0 Axes>
Image

Figure 3: Above, the calexp image with the subtracted sky background added back in. Note that the scale bar now ranges from 3500 to 4200 counts.

Clean up. Do not use the calexp + bkgd in the sections below. A fresh version of the calexp is retrieved in Section 3.1.2.

In [13]:
del calexp, bkgd

3. Source detection, deblending, and measurement¶

The LSST Science Pipelines' source detection, deblending, and measurement tasks were imported in the first code cell.

More information can be found at pipelines.lsst.io, either by scrolling down and browsing the available modules or by using the search bar at upper left.

The schema describes the output properties that will be measured for each source. The schema needs to be passed to all of the tasks, each of which will add columns to it when it runs.

Create a minimal schema for a source table.

In [14]:
schema = afwTable.SourceTable.makeMinimalSchema()
print(schema)

Schema(
    (Field['L'](name="id", doc="unique ID"), Key<L>(offset=0, nElements=1)),
    (Field['Angle'](name="coord_ra", doc="position in ra/dec"), Key<Angle>(offset=8, nElements=1)),
    (Field['Angle'](name="coord_dec", doc="position in ra/dec"), Key<Angle>(offset=16, nElements=1)),
    (Field['L'](name="parent", doc="unique ID of parent source"), Key<L>(offset=24, nElements=1)),
)

Create a container which will be used to record metadata about algorithm execution.

In [15]:
algMetadata = dafBase.PropertyList()

3.1. Configure the tasks¶

Each task possesses an associated configuration class. The properties of these configuration classes can be determined from the classes themselves.

Option: uncomment the following line to view the help documentation for the CharacterizeImageTask configuration. Replace CharacterizeImageTask with another task name to view its help documentation.

In [16]:
# help(CharacterizeImageTask.ConfigClass())
In [17]:
# help(SourceDetectionTask.ConfigClass())

Set the basic configuration parameters and instantiate each task in turn.

CharacterizeImageTask characterizes the image properties (e.g., the PSF).

Set the configuration to only do one PSF iteration.

In [18]:
config = CharacterizeImageTask.ConfigClass()
config.psfIterations = 1
charImageTask = CharacterizeImageTask(config=config)
del config

SourceDetectionTask detects sources.

Set the configuration for the detection threshold to 4 (the default for LSST catalogs is 5) in order to detect sub-threshold (faint) sources. Set the detection threshold type to be stdev, which means the threshold is applied to the image's standard deviation. In other words, this configuration will detect sources at $4\sigma$.

In [19]:
config = SourceDetectionTask.ConfigClass()
config.thresholdValue = 4
config.thresholdType = "stdev"
sourceDetectionTask = SourceDetectionTask(schema=schema, config=config)
del config

SourceDeblendTask deblends overlapping sources into "parent" and "child" sources. No configuration is needed for this tutorial.

In [20]:
sourceDeblendTask = SourceDeblendTask(schema=schema)

SingleFrameMeasurementTask measures the properties of the deblended sources. No special configuration is needed for this tutorial.

In [21]:
config = SingleFrameMeasurementTask.ConfigClass()
sourceMeasurementTask = SingleFrameMeasurementTask(schema=schema,
                                                   config=config,
                                                   algMetadata=algMetadata)
del config

Create a SourceTable called tab using the minimal schema. It will hold the output of the analysis in the following section.

Then, as schema is no longer needed and a fresh minimal schema is made in Section 5, delete this one.

In [22]:
tab = afwTable.SourceTable.make(schema)
del schema

3.1.1. Re-configure tasks and explore methods¶

The configuration parameters cannot be changed after the task is constructed. To change a configuration parameter, redefine it and then also redefine the task.

For example, to change the number of PSF iterations done during image characterization:

config.psfIterations = 3
charImageTask = CharacterizeImageTask(config=config)

Option: To explore a task, use the help function on the task or any of its methods.

In [23]:
# help(charImageTask)
In [24]:
# help(charImageTask.writeSchemas)

Option: To see a pop-up list of all methods for a given task, uncomment this line, place the cursor after the period, and press tab. Then re-comment the line. Do not execute the cell.

In [25]:
# charImageTask.

3.1.2. Retrieve a fresh version of the image¶

Retrieve a fresh, unaltered calexp image from the butler, and again clear the "DETECTED" information from the mask plane.

In [26]:
calexp = butler.get('calexp', dataId=dataId)
In [27]:
calexp.mask.removeAndClearMaskPlane('DETECTED')

3.2. Image characterization¶

Characterize the calexp image. This calculates various global properties of the image, such as the PSF FWHM (full-width half-max).

In [28]:
result = charImageTask.run(calexp)

lsst.characterizeImage INFO: PSF estimation initialized with 'simple' PSF
lsst.characterizeImage.repair INFO: Identified and kept 0 cosmic rays.
lsst.characterizeImage.detection INFO: Setting factor for negative detections equal to that for positive detections: 1.000000
lsst.characterizeImage.detection INFO: Detected 249 positive peaks in 212 footprints and 0 negative peaks in 0 footprints to 50 +ve and 50 -ve sigma
lsst.characterizeImage.detection INFO: Resubtracting the background after object detection
lsst.characterizeImage.measurement INFO: Measuring 212 sources (212 parents, 0 children) 
lsst.characterizeImage.measurePsf INFO: Measuring PSF
lsst.characterizeImage.measurePsf INFO: PSF star selector found 124 candidates
lsst.characterizeImage.measurePsf.reserve INFO: Reserved 0/124 sources
lsst.characterizeImage.measurePsf INFO: Sending 124 candidates to PSF determiner
lsst.characterizeImage.measurePsf INFO: PSF determination using 123/124 stars.
lsst.characterizeImage INFO: iter 1; PSF sigma=2.0684, dimensions=(41, 41); median background=0.83
lsst.characterizeImage.repair INFO: Identified and interpolated over 0 cosmic rays.
lsst.characterizeImage.maskStreaks INFO: The Kernel Hough Transform detected 0 line(s)
lsst.characterizeImage.measurement INFO: Measuring 212 sources (212 parents, 0 children) 
lsst.characterizeImage.measureApCorr INFO: Measuring aperture corrections for 2 flux fields
lsst.characterizeImage.measureApCorr.sourceSelector INFO: Selected 64/212 sources
lsst.characterizeImage.measureApCorr INFO: Aperture correction for base_GaussianFlux from 63 stars: MAD 0.002821, RMS 0.003946
lsst.characterizeImage.measureApCorr INFO: Aperture correction for base_PsfFlux from 63 stars: MAD 0.003641, RMS 0.004337
lsst.characterizeImage.applyApCorr INFO: Applying aperture corrections to 2 instFlux fields

The result is a structure, lsst.pipe.base.struct.Struct.

A Struct is just a generalized container for storing multiple output components and accessing them as attributes.

In [29]:
type(result)
Out[29]:
lsst.pipe.base.struct.Struct

The content of this Struct can be investigated with the getDict method.

In [30]:
for k, v in result.getDict().items():
    print(k, type(v))

exposure <class 'lsst.afw.image._exposure.ExposureF'>
sourceCat <class 'lsst.afw.table.SourceCatalog'>
background <class 'lsst.afw.math._backgroundList.BackgroundList'>
psfCellSet <class 'NoneType'>
characterized <class 'lsst.afw.image._exposure.ExposureF'>
backgroundModel <class 'lsst.afw.math._backgroundList.BackgroundList'>

Option: the Struct can be displayed but this is not very useful.

In [31]:
# result

As an example of how to interact with the results of image characterization, get the PSF at a specific pixel location and print the FWHM. The factor of 2.355 converts from the standard deviation to the FWHM.

In [32]:
x, y = 1700, 2100
point = geom.Point2D(x, y)
psf = calexp.getPsf()
sigma = psf.computeShape(point).getDeterminantRadius()
pixelScale = calexp.getWcs().getPixelScale().asArcseconds()
print('PSF FWHM = {:.2f} arcsec'.format(sigma*pixelScale*2.355))

del x, y, point, psf, sigma, pixelScale

PSF FWHM = 0.97 arcsec

3.3. Source detection¶

Run source detection.

In [33]:
result = sourceDetectionTask.run(tab, calexp)

lsst.sourceDetection INFO: Setting factor for negative detections equal to that for positive detections: 1.000000
lsst.sourceDetection INFO: Detected 2292 positive peaks in 1983 footprints and 109 negative peaks in 108 footprints to 4 +ve and 4 -ve sigma
lsst.sourceDetection INFO: Resubtracting the background after object detection

Print the number of positive peaks detected.

In [34]:
result.numPosPeaks
Out[34]:
2292

With thresholdValue set to 4, over 2000 sources are detected.

Had the thresholdValue been set instead to 10, there would have been about 940 peaks detected.

Create sources to hold the detected sources from the results.

In [35]:
sources = result.sources

Option: Display sources to see that most of the columns are NaN at first.

In [36]:
# sources

Option: Write the sources to a FITS table, or save the calexp as a FITS file, in the home directory.

In [37]:
# home_directory = '/home/' + getpass.getuser() + '/'
# sources.writeFits(home_directory + "DP02_NB05_sources.fits")
# calexp.writeFits(home_directory + "DP02_NB05_calexp.fits")

3.4. Source debelending and measurement¶

Run the SourceDeblendTask and SingleFrameMeasurementTask.

In [38]:
sourceDeblendTask.run(calexp, sources)

lsst.sourceDeblend INFO: Deblending 1983 sources
lsst.sourceDeblend INFO: Deblended: of 1983 sources, 227 were deblended, creating 520 children, total 2503 sources
In [39]:
sourceMeasurementTask.run(measCat=sources, exposure=calexp)

lsst.measurement INFO: Measuring 2503 sources (1983 parents, 520 children)

Ensure that the sources catalog is contiguous (sequential) in memory by making a copy.

In [40]:
sources = sources.copy(True)

Option: Display the sources catalog as an astropy table, which has a nice, human-readable output format.

In [41]:
# sources.asAstropy()

3.5. Visualize sources on an image cutout¶

Overplot the detected sources on a cutout image of the calexp using afwDisplay.

There are several ways to create image cutouts. This tutorial demonstrates the use of the Factory method for a calexp.

The arguments that must be passed to the Factory method are:

calexp : the ExposureF to make the cutout from
bbox   : the desired bounding box of the cutout
origin : the image pixel origin, local to the cutout array
deep   : if True, copy the data rather than passing by reference

Define the bbox region for a cutout that is centered on pixel coordinate 1700, 2100 and is 400 pixels wide and 400 pixels high.

In [42]:
x, y = 1700, 2100
w, h = 400, 400
xmin, ymin = (x - w)//2, (y - h)//2
bbox = geom.Box2I()
bbox.include(geom.Point2I(xmin, ymin))
bbox.include(geom.Point2I(xmin + w, ymin + h))
del x, y, w, h, xmin, ymin

Option: Display the bbox information.

In [43]:
# bbox

An alternative statement to create bbox is:

bbox = geom.Box2I(geom.Point2I(xmin, ymin), geom.Extent2I(w, h))

Generate the cutout.

In [44]:
cutout = calexp.Factory(calexp, bbox, origin=afwImage.LOCAL, deep=False)

Display the cutout and overplot sources as orange circles.

Use display buffering while overplotting sources to avoid re-drawing the image after each source is plotted.

In [45]:
fig = plt.figure()
afw_display = afwDisplay.Display()
afw_display.scale('asinh', 'zscale')
afw_display.mtv(cutout.image)
plt.gca().axis('off')
with afw_display.Buffering():
    for s in sources:
        afw_display.dot('o', s.getX(), s.getY(), size=20, ctype='orange')
plt.show()

<Figure size 800x800 with 0 Axes>
Image

Figure 4: Above, a cutout of the calexp with about 20 detected sources circled in orange. Since a low detection threshold was used, some orange circles do not have a clearly visible source in them.

Clean up. Keep sources and result for the next sections.

In [46]:
del bbox, cutout, tab

4. Footprints¶

Object footprints are an integral component of the high-level CCD processing tasks (e.g., detection, measurement, and deblending).

To quote Bosch et al. (2017): "Footprints record the exact above-threshold detection region on a CCD. These are similar to SExtractor’s “segmentation map", in that they identify which pixels belong to which detected objects."

This quote draws an analogy between footprints and segmentation maps because they both identify pixels with values above some threshold.

The result of the SourceDetectionTask stores the footprints associated with detected objects.

Grab the above-threshold footprints that were detected, and assign them to the variable fps.

In [47]:
fpset = result.positive
fps = fpset.getFootprints()

Get a rough view of the first source's footprint from its span.

In [48]:
fps[0].getSpans()
Out[48]:
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Above, the footprint is represented by the 1s in the array.

Note that the first row of the span array will be the bottom row of the image, and so is "upside down" compared to the heavy footprint visualization below.

4.1. Heavy Footprints¶

The footprint span above indicates which pixels are included in the footprint (the 1s).

A heavy fooprint includes the pixel value data from the image.

Demonstrate that the fps footprint is not heavy.

In [49]:
fps[0].isHeavy()
Out[49]:
False

Make all of the footprints heavy and define hfps to hold the heavy footprints.

In [50]:
fpset.makeHeavy(calexp.getMaskedImage())
hfps = fpset.getFootprints()

The getImageArray method returns a flattened, 1-dimensional array of pxiels from the footprint.

Option: Uncomment the line below to print the array and see that it contains pixel values.

In [51]:
# hfps[0].getImageArray()

Use the span's method to unflatten the image array, and display the heavy footprint (matplotlib colormap options).

In [52]:
fig = plt.figure()
plt.imshow(fps[0].getSpans().unflatten(hfps[0].getImageArray()),
           cmap='bone', origin='lower')
plt.show()
Image

The heavy footprint also comes with a 1-dimensional mask array.

Option: Display the 1-d mask array.

In [53]:
# hfps[0].getMaskArray()

Review the dictionary for the mask plane, which maps mask value to the reason the pixel is masked.

In [54]:
calexp.getMask().getMaskPlaneDict()
Out[54]:
{'BAD': 0,
 'CR': 3,
 'CROSSTALK': 9,
 'DETECTED': 5,
 'DETECTED_NEGATIVE': 6,
 'EDGE': 4,
 'INTRP': 2,
 'NOT_DEBLENDED': 10,
 'NO_DATA': 8,
 'SAT': 1,
 'STREAK': 12,
 'SUSPECT': 7,
 'UNMASKEDNAN': 11}

The values are the exponent of the bitmask.

Pixels only marked detected will be 2^5 = 32.

Pixels that are both on the edge of the original image and detected will be 2^5 + 2^4 = 48.

Visualize mask array.

To show the colorbar, create a new axis cax on the right side of the image display, with a width of 5% and a padding between display and colorbar of 1%. Add ticks to the colorbar at values of 32 and 48, to correspond to the two mask values of interest ("detected" and "detected + edge").

In [55]:
plt.figure()
im = plt.imshow(fps[0].getSpans().unflatten(hfps[0].getMaskArray()),
                origin='lower')
ax = plt.gca()
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad="1%")
plt.colorbar(im, cax=cax, ticks=[0, 32, 32+16])
plt.show()
Image

Figure 5: Above, the footprint's mask array reveals which pixels are considered (green) but which are also at the edge of the calexp (yellow).

Clean up.

In [56]:
del fpset, fps, hfps

5. Forced photometry¶

The LSST Science Pipelines have a task for forced photometry, called ForcedMeasurementTask.

This task will use the PSF of the image to measure flux at any user-supplied coordinates. As it is not a fit, a flux is returned whether an astronomical source is there or not. Forced photometry can be very useful for including sources below the detection threshold (or non-detections) in analyses.

5.1. Configure the task¶

To start, create a minimal schema.

In [57]:
schema = afwTable.SourceTable.makeMinimalSchema()

At this time, the ForcedMeasurementTask expects a few additional columns related to the centroids and shapes to pre-exist, so add them to the minimal schema. These might change in the future.

In [58]:
alias = schema.getAliasMap() 
x_key = schema.addField("centroid_x", type="D")
y_key = schema.addField("centroid_y", type="D")
alias.set("slot_Centroid", "centroid")

xx_key = schema.addField("shape_xx", type="D")
yy_key = schema.addField("shape_yy", type="D")
xy_key = schema.addField("shape_xy", type="D")
alias.set("slot_Shape", "shape")

Add a type flag, which will be 0 for sub-threshold (faint) sources and 1 for random "off-source" coordinates.

In [59]:
type_key = schema.addField("type_flag", type="F")

Create an empty forcedSource catalog using the schema. It will have 0 rows, to start.

In [60]:
forcedSource = afwTable.SourceCatalog(schema)
print(len(forcedSource))

0

Configure the ForcedMeasurementTask. At this time, few plugins beyond the default are required, as is the doReplaceWithNoise configuration. These might change in the future.

In [61]:
config = ForcedMeasurementTask.ConfigClass()
config.copyColumns = {}
config.plugins.names = [
    "base_PixelFlags",
    "base_TransformedCentroid",
    "base_PsfFlux",
    "base_TransformedShape"
]
config.doReplaceWithNoise = False

Define the task using the schema and configuration.

In [62]:
forcedMeasurementTask = ForcedMeasurementTask(schema, config=config)
del config

5.2. Define the coordinates¶

Scientific motivations for forced photometry include measurement of sub-threshold (faint) sources or a need for an upper-limit on the flux of a time-varying source.

The potential drawback of forced photometry is that all the flux within the area of the PSF is included, which means nearby sources can contaminate (add excess flux to) the forced photometry. This will be demonstrated in Section 5.4.

For this example, create a list of sub-threshold (faint) coordinates and random "off-source" coordinates which are not near any detected sources. Use a region within 2 arcmin of the center of the field.

5.2.1. Sub-threshold (faint) sources¶

Create numpy arras of the right ascension and declination, x- and y-pixel coordinates, and instrumental flux and its error from the list of detected sources. Calculate the signal-to-noise ratio.

In [63]:
srcs_ra = np.asarray(np.rad2deg(sources['coord_ra']), dtype='float')
srcs_dec = np.asarray(np.rad2deg(sources['coord_dec']), dtype='float')
srcs_xpix = np.asarray(sources['base_SdssCentroid_x'], dtype='float')
srcs_ypix = np.asarray(sources['base_SdssCentroid_y'], dtype='float')
srcs_flux = np.asarray(sources['base_PsfFlux_instFlux'], dtype='float')
srcs_fluxe = np.asarray(sources['base_PsfFlux_instFluxErr'], dtype='float')
srcs_snr = srcs_flux/srcs_fluxe

Use the bbox and wcs to get the central coordinates of the calexp, convert from pixels to sky coordinates, and calculate offsets (s_dra and s_ddec) between the sources and the image center.

In [64]:
bbox = calexp.getBBox()
wcs = calexp.getWcs()
center = wcs.pixelToSky(bbox.centerX, bbox.centerX)
cra = center.getRa().asDegrees()
cdec = center.getDec().asDegrees()
srcs_dra = np.abs(srcs_ra - cra)
srcs_ddec = np.abs(srcs_dec - cdec)

Define the maximum offset from the image center to be 2 arcminutes (in units of degrees).

Find all sources with 4 $<$ SNR $<$ 5 (i.e., sub-threshold sources) within the central region.

Define new arrays to hold only the information of the sources to be used for forced photometry, and clean up.

In [65]:
dmax = 2.0 / 60.0
tx = np.where((srcs_snr > 4) & (srcs_snr < 5)
              & (srcs_dra < dmax) & (srcs_ddec < dmax))[0]
print(len(tx))

s_ra = srcs_ra[tx]
s_dec = srcs_dec[tx]
s_xpix = srcs_xpix[tx]
s_ypix = srcs_ypix[tx]
s_flux = srcs_flux[tx]

del tx, dmax
del srcs_ra, srcs_dec, srcs_xpix, srcs_ypix, srcs_flux, srcs_fluxe, srcs_snr
del srcs_dra, srcs_ddec

32

Add the coordinates of these 32 sub-threshold sources to the forcedSource catalog.

Sky coordinates must be in radians.

In [66]:
for i in range(len(s_ra)):
    sourceRec = forcedSource.addNew()
    coord = geom.SpherePoint(np.radians(s_ra[i]) * geom.radians,
                             np.radians(s_dec[i]) * geom.radians)
    sourceRec.setCoord(coord)
    sourceRec[x_key] = s_xpix[i]
    sourceRec[y_key] = s_ypix[i]
    sourceRec[type_key] = 0
    del sourceRec, coord

print(len(forcedSource))

32

5.2.2. Random "off-source" coordinates¶

Create 32 random coordinates in the central region of the image. Ensure they are offset by at least 20 pixels from the nearest detected source.

In [67]:
pixelScale = calexp.getWcs().getPixelScale().asArcseconds()
radius_pix = int(2.0 * 60.0 / pixelScale)
center_x = bbox.centerX
center_y = bbox.centerY

tempx = []
tempy = []
while len(tempx) < 32:
    r = np.random.choice(radius_pix)
    theta = np.random.rand() * 2 * np.pi
    x = center_x + r * np.cos(theta)
    y = center_y + r * np.sin(theta)
    off = np.min(np.sqrt((x - s_xpix)**2 + (y - s_ypix)**2))
    if off > 20:
        tempx.append(x)
        tempy.append(y)
    del r, theta, x, y, off

r_xpix = np.asarray(tempx)
r_ypix = np.asarray(tempy)
del pixelScale, radius_pix, center_x, center_y, tempx, tempy

Use the image's WCS to convert the random pixels to sky coordinates.

In [68]:
r_ra = np.zeros(len(r_xpix), dtype='float')
r_dec = np.zeros(len(r_ypix), dtype='float')
for i in range(len(r_xpix)):
    r_skycoord = wcs.pixelToSky(r_xpix[i], r_ypix[i])
    r_ra[i] = r_skycoord.getRa().asDegrees()
    r_dec[i] = r_skycoord.getDec().asDegrees()
    del r_skycoord

Add the random "off-source" coordinates to the forcedSource catalog.

In [69]:
for i in range(len(r_xpix)):
    sourceRec = forcedSource.addNew()
    coord = geom.SpherePoint(np.radians(r_ra[i]) * geom.radians,
                             np.radians(r_dec[i]) * geom.radians)
    sourceRec.setCoord(coord)
    sourceRec[x_key] = r_xpix[i]
    sourceRec[y_key] = r_ypix[i]
    sourceRec[type_key] = 1
    del sourceRec, coord

print(len(forcedSource))

64

Option: View the forcedSource table.

In [70]:
# forcedSource

5.3. Run the task¶

Create an empty catalog to store forced measurements.

In [71]:
forcedMeasCat = forcedMeasurementTask.generateMeasCat(calexp, forcedSource,
                                                      calexp.getWcs())

Option to review the forcedMeasCat. Before the task is run, all the columns are empty (or default values).

In [72]:
# forcedMeasCat

Run the forced measurements.

In [73]:
forcedMeasurementTask.run(forcedMeasCat, calexp, forcedSource, calexp.getWcs())

lsst.measurement INFO: Performing forced measurement on 64 sources

5.4. Explore the results¶

Get the results as an astropy table.

In [74]:
table = forcedMeasCat.asAstropy()

Add the type_flag column from forcedSource to table.

In [75]:
temp = forcedSource.asAstropy()
table['type_flag'] = temp['type_flag']
del temp

Option: View the filled forcedMeasCat results as the astropy table.

In [76]:
# table

Overplot the coordinates used for forced photometry on the image.

In [77]:
plt.figure(figsize=(8, 8))
afw_display = afwDisplay.Display()
afw_display.scale('linear', 'zscale')
afw_display.mtv(calexp.image)
plt.gca().axis('off')
with afw_display.Buffering():
    for i in range(len(table)):
        clr = 'orange'
        if table['type_flag'][i] == 1:
            clr = 'cyan'
        afw_display.dot('o', table['slot_Centroid_x'][i],
                        table['slot_Centroid_y'][i],
                        size=20, ctype=clr)
plt.xlim([np.min(s_xpix)-100, np.max(s_xpix)+100])
plt.ylim([np.min(s_ypix)-100, np.max(s_ypix)+100])
plt.show()

<Figure size 800x800 with 0 Axes>
Image

Figure 6: Above, a zoom-in to the central region of the image with sub-threshold (faint) source coodinates (orange) and random "off-source" coordinates (cyan) showing the locations where forced photometry was done.

Compare the forced photometry fluxes with the measured fluxes for sub-threshold (faint) sources.

In [78]:
fig = plt.figure(figsize=(6, 4))
plt.plot([2500, 3000], [2500, 3000], lw=0.5, color='black')
for i in range(32):
    plt.plot(s_flux[i], table['base_PsfFlux_instFlux'][i],
             'o', color='darkorange', mew=0, alpha=0.5)
plt.xlabel('detection flux')
plt.ylabel('forced PSF flux')
plt.show()
Image

Figure 7: Above, the original detection instrumental flux (in counts) versus the forced PSF flux measurements for the 32 sources which were originally detected with 4 $<$ SNR $<$ 5. For most sources the fluxes agree, but in some cases the forced PSF flux is greater.

Cases in which the forced PSF flux is greater is likely due to contamination from a nearby source. This contamination is mitigated when source deblending and measurement are done together, but forced photometry includes all the flux in the area of the PSF.

Again overplot sources on the image, but color as yellow the ones with a flux differential greater than 20 counts.

In [79]:
plt.figure(figsize=(8, 8))
afw_display = afwDisplay.Display()
afw_display.scale('linear', 'zscale')
afw_display.mtv(calexp.image)
plt.gca().axis('off')
with afw_display.Buffering():
    for i in range(len(table)):
        clr = 'orange'
        if table['type_flag'][i] == 0:
            if table['base_PsfFlux_instFlux'][i] > s_flux[i] + 20:
                clr = 'yellow'
            afw_display.dot('o', table['slot_Centroid_x'][i],
                            table['slot_Centroid_y'][i],
                            size=20, ctype=clr)
plt.xlim([np.min(s_xpix)-100, np.max(s_xpix)+100])
plt.ylim([np.min(s_ypix)-100, np.max(s_ypix)+100])
plt.show()

<Figure size 800x800 with 0 Axes>
Image

Figure 8: Above, the yellow circles mark locations of sub-threshold (faint) sources for which the forced photometry had a flux excess, and for at least five the contaminating source is visible.

As a final exploration, plot the histograms of the forced PSF fluxes for the sub-threshold (faint) sources and the random "off-source" coordinates.

In [80]:
fig = plt.figure(figsize=(6, 4))
tx1 = np.where(table['type_flag'] == 0)[0]
tx2 = np.where(table['type_flag'] == 1)[0]
plt.hist(table['base_PsfFlux_instFlux'][tx1], bins=10, histtype='step',
         color='darkorange', label='sub-threshold source')
plt.hist(table['base_PsfFlux_instFlux'][tx2], bins=40, histtype='step',
         color='dodgerblue', label='random coordinate')
plt.xlabel('forced PSF flux')
plt.ylabel('number of forced sources')
plt.legend(loc='upper left')
plt.show()
Image